#-----------------------------------------------------------------------------------------------------------------
#-------------------------------------------Main Tables and Figures-----------------------------------------------
#-----------------------------------------------------------------------------------------------------------------

#R version 4.0.3 (2020-10-10)

## Required packages (May need to restart R program after installation)
install.packages("plyr")
install.packages("dplyr")
install.packages("estimatr")
install.packages("multcomp") # may require "mvtnorm", "survival", "TH.data", "MASS" packages
install.packages("readstata13")
install.packages("ggplot2")


## Load libraries
library(readstata13) # for loading data
library(plyr) # for data processing
library(dplyr) # for data processing
library(estimatr) # for estimation
library(multcomp) # for general linear hypotheses testing
library(ggplot2) # for plot


rm(list = ls())

## Set working directory
#Please set the working directory to the folder where the readme.txt file is located
setwd("xxx/xxx/xxx")

## load college survey data
survey2nd <-read.dta13("data/collegesurvey.dta")

## generate factor treatment variable
survey2nd$groupc_fct <- factor(survey2nd$qc_group, labels = c ("Control", "Social", "Political", "Social.Political"))


##-----------------Table 1. Experimental Design for Attitude toward the SCS---------------------
## Count observations for Table 1
table(survey2nd$groupc_fct)
#         Control           Social        Political Social.Political 
#             204              164              198              181 


## create a file to standardize variables
survey2nd_std <- survey2nd

## generate a social-order treatment dummy for factorial regression
survey2nd_std <- survey2nd_std %>% 
  mutate(treatcvorg = treat_cv) %>%
  as.data.frame()

survey2nd_std <- survey2nd_std %>%
  mutate(treat_cv = replace(treat_cv, treat_cv_po==1, 1) ) %>%
  as.data.frame()

#table(survey2nd_std$treat_cv)
#table(survey2nd_std$treatcvorg)

## generate a repression treatment dummy for factorial regression
survey2nd_std <- survey2nd_std %>% 
  mutate(treatpoorg = treat_po) %>%
  as.data.frame()

survey2nd_std <- survey2nd_std %>%
  mutate(treat_po = replace(treat_po, treat_cv_po==1, 1) ) %>%
  as.data.frame()

#table(survey2nd_std$treat_po)
#table(survey2nd_std$treatpoorg)


## standardize variables
survey2nd_std$scsapproval <- scale(survey2nd_std$sc_legit)
survey2nd_std$distrust <- scale(survey2nd_std$op_takeother)
survey2nd_std$distrustexp <- scale(survey2nd_std$op_takeself)
survey2nd_std$propaganda <- scale(survey2nd_std$scinfo_state)
survey2nd_std$femalestd <- scale(survey2nd_std$female)
survey2nd_std$agestd <- scale(survey2nd_std$age)
survey2nd_std$incomestd <- scale(survey2nd_std$faminc)
survey2nd_std$ccpstd <- scale(survey2nd_std$party_status)
survey2nd_std$mediastd <- scale(survey2nd_std$media_phone)
survey2nd_std$behaviorstd <- scale(survey2nd_std$behavior)
survey2nd_std$servicestd <- scale(survey2nd_std$par_org)
survey2nd_std$scoreselfstd <- scale(survey2nd_std$sc_selfscore)

survey2nd_std$treatcvstd <- scale(survey2nd_std$treat_cv)
survey2nd_std$treatpostd <- scale(survey2nd_std$treat_po)
survey2nd_std$treatcvpostd <- scale(survey2nd_std$treat_cv_po)



## generate a variable indicating state media as the only source of information about the SCS
survey2nd_std$infostateonly <- survey2nd_std$scinfo_state

survey2nd_std <- survey2nd_std %>%
  mutate(infostateonly = replace(infostateonly, scinfo_nonstate==1, 0) ) %>%
  as.data.frame()

# check the variable
table(survey2nd_std$infostateonly)
#  0   1 
#557 180 

## main model (full sample)
mod.exp <- lm_robust(scsapproval ~ treatcvstd + treatpostd + treatcvpostd +  factor(univ),
                     data = survey2nd_std, se_type = "stata")
mod.exp$nobs
#[1] 737

## model: nonstate media as sources
mod.expinfo <- lm_robust(scsapproval ~ treatcvstd + treatpostd + treatcvpostd + factor(univ), 
                         data = subset(survey2nd_std, infostateonly == 0), se_type = "stata")
mod.expinfo$nobs
#[1] 557

## model: state media as the only source
mod.expnoinfo <- lm_robust(scsapproval ~ treatcvstd + treatpostd + treatcvpostd + factor(univ), 
                           data = subset(survey2nd_std, infostateonly == 1), se_type = "stata")
mod.expnoinfo$nobs
#[1] 180

## test linear combinations for interaction models
testpo <- summary(glht(mod.exp, linfct = c("treatpostd+treatcvpostd = 0")))
testcv <- summary(glht(mod.exp, linfct = c("treatcvstd+treatcvpostd = 0")))

testpoi <- summary(glht(mod.expinfo, linfct = c("treatpostd+treatcvpostd = 0")))
testcvi <- summary(glht(mod.expinfo, linfct = c("treatcvstd+treatcvpostd = 0")))

testponi <- summary(glht(mod.expnoinfo, linfct = c("treatpostd+treatcvpostd = 0")))
testcvni <- summary(glht(mod.expnoinfo, linfct = c("treatcvstd+treatcvpostd = 0")))

## extract test statistics and make tables for plots
tab.exp0 <- rbind(
  cbind(testpo$test$coefficients, testpo$test$coefficients - 1.96* testpo$test$sigma, testpo$test$coefficients + 1.96* testpo$test$sigma),
  cbind(testcv$test$coefficients, testcv$test$coefficients - 1.96* testcv$test$sigma, testcv$test$coefficients + 1.96* testcv$test$sigma)
)

tab.expinfo0 <- rbind(
  cbind(testpoi$test$coefficients, testpoi$test$coefficients - 1.96* testpoi$test$sigma, testpoi$test$coefficients + 1.96* testpoi$test$sigma),
  cbind(testcvi$test$coefficients, testcvi$test$coefficients - 1.96* testcvi$test$sigma, testcvi$test$coefficients + 1.96* testcvi$test$sigma)
)

tab.expnoinfo0 <- rbind(
  cbind(testponi$test$coefficients, testponi$test$coefficients - 1.96* testponi$test$sigma, testponi$test$coefficients + 1.96* testponi$test$sigma),
  cbind(testcvni$test$coefficients, testcvni$test$coefficients - 1.96* testcvni$test$sigma, testcvni$test$coefficients + 1.96* testcvni$test$sigma)
)

tab.exp <- cbind(mod.exp$coefficients[4:2],mod.exp$conf.low[4:2],mod.exp$conf.high[4:2])
tab.expinfo <- cbind(mod.expinfo$coefficients[4:2],mod.expinfo$conf.low[4:2],mod.expinfo$conf.high[4:2])
tab.expnoinfo <- cbind(mod.expnoinfo$coefficients[4:2],mod.expnoinfo$conf.low[4:2],mod.expnoinfo$conf.high[4:2])

nas<- cbind(NA,NA,NA)

tab.exp <-rbind(tab.exp0, nas, tab.exp)
tab.expinfo <-rbind(tab.expinfo0, nas, tab.expinfo)
tab.expnoinfo <-rbind(tab.expnoinfo0, nas, tab.expnoinfo)

## function for word wrap in figures
wordwrap<-function(x,len) paste(strwrap(x,width=len),collapse="\n") 

#----------------------Figure 3. Information Treatment Effects: Full Sample-------------------
pdf("outputmain/expttl.pdf", width=6.5, height=5) 
par(mar=c(3.9, 9, 1, 1) + 1)
plot(tab.exp[,1],1:6,xlim=c(-0.6,0.6),ylim=c(0.5, 6.5),main = "",ylab = "",cex.lab=1,cex.axis=1, yaxt="n",xlab="Coefficients of Standardized Variables with 95% CIs")
points(tab.exp[1:2,1],1:2, pch=1)
points(tab.exp[4:6,1],4:6, pch=16, cex=1.5)
axis(2,at=c(6,5,4,2,1),labels=sapply(c("Social Order", "Political Repression", "Order & Repression", "Order+Interaction", "Repression+Interaction"
),wordwrap, len=20),cex.axis=1,mgp=c(2,1,0),las=1)
#mtext("Treatment",side=1,line=6.5,cex=1.5)
abline(v=0,lty=2,col="red")
abline(h=3)
legend(-0.655, 6.65, legend="Point Estimates:", box.lty=0)
legend(-0.655, 2.85, legend="Linear Combinaions:", box.lty=0)
for(i in 1:2) segments(tab.exp[i,2],i,tab.exp[i,3],i,lwd=1)
for(i in 4:6) segments(tab.exp[i,2],i,tab.exp[i,3],i,lwd=1.5)
dev.off()

#----------------------Figure 4. Information Treatment Effects, By Information Sources-------------------
pdf("outputmain/expttlcomp.pdf", width=6.5, height=5) 
par(mar=c(3.9, 9, 1, 1) + 1)
plot(tab.expinfo[,1],1:6,xlim=c(-0.6,0.6),ylim=c(0.41, 6.41),main = "",ylab = "",pch="",cex.lab=1,cex.axis=1, yaxt="n",xlab="Coefficients of Standardized Variables with 95% CIs")
points(tab.expinfo[1:2,1],1:2, pch=1)
points(tab.expinfo[4:6,1],4:6, pch=16, cex=1.5)
axis(2,at=c(6,5,4,2,1),labels=sapply(c("Social Order", "Political Repression", "Order & Repression", "Order+Interaction", "Repression+Interaction"
),wordwrap, len=20),cex.axis=1,mgp=c(2,1,0),las=1)
#mtext("Treatment",side=1,line=6.5,cex=1.5)
abline(v=0,lty=2,col="red")
abline(h=3)
legend(-0.655, 6.65, legend="Point Estimates:", box.lty=0)
legend(-0.655, 2.85, legend="Linear Combinaions:", box.lty=0)
for(i in 1:2) segments(tab.expinfo[i,2],i,tab.expinfo[i,3],i,lwd=1)
for(i in 4:6) segments(tab.expinfo[i,2],i,tab.expinfo[i,3],i,lwd=1.5)
par(new=TRUE)
plot(tab.expnoinfo[,1],1:6,xlim=c(-0.6,0.6),ylim=c(0.59, 6.59),main = "",ylab = "",pch="",cex.lab=1,cex.axis=1, yaxt="n",xlab="Coefficients of Standardized Variables with 95% CIs")
points(tab.expnoinfo[1:2,1],1:2, pch=2)
points(tab.expnoinfo[4:6,1],4:6, pch=17, cex=1.5)
for(i in 1:2) segments(tab.expnoinfo[i,2],i,tab.expnoinfo[i,3],i,lwd=1)
for(i in 4:6) segments(tab.expnoinfo[i,2],i,tab.expnoinfo[i,3],i,lwd=1.5)
legend(0.15, 1.5, legend=c("More Informed", "Less Informed"), pch=1:2, lty=1, lwd=0.75, cex=0.75,pt.cex=0.75, box.lty=0)
legend(0.15, 4.1, legend=c("More Informed", "Less Informed"), pch=16:17, lty=1, lwd=1, cex=0.75,pt.cex=1, box.lty=0)
dev.off()




## drop insincere respondents for observational studies. Note: the results are very similar when keeping these respondents
nrow(survey2nd_std)
#[1] 754
survey2ndobs_std <- subset(survey2nd_std, provcppcc < 1) # drop respondents who failed the attention check question.
nrow(survey2ndobs_std)
#[1] 738


## observational study model
mod.approvalexp2 <- lm_robust(scsapproval ~ propaganda + distrust + distrustexp +
                                agestd + femalestd + incomestd + ccpstd + behaviorstd +servicestd + factor(univ), cluster = univ,
                              data = survey2ndobs_std, se_type = "stata")
mod.approvalexp2$nobs
#[1] 665
tab.approvalexp2 <- cbind(mod.approvalexp2$coefficients[10:2],mod.approvalexp2$conf.low[10:2],mod.approvalexp2$conf.high[10:2])


#----------------------Figure 5. Sources of Support for SCSs: College Field Survey-------------------
pdf("outputmain/approvalexp2.pdf", width=7.9, height=5) 
par(mar=c(3.9, 15, 1, 1) + 1)
plot(tab.approvalexp2[,1],1:9,xlim=c(-0.5,0.5),ylim=c(0.5, 9.5),main = "",ylab = "",pch="",cex.lab=1,cex.axis=1, yaxt="n",xlab="Coefficients of Standardized Variables with 95% CIs")
points(tab.approvalexp2[1:8,1],1:8, pch=1)
points(tab.approvalexp2[9,1],9, pch=16, cex=1.5)
axis(2,9:1,labels=sapply(c("SCS Info. Source: State Media", "Perceived Distrust in Society", "Personal Distrust", "Age", "Female", "Income", "CCP Member", "Social Conformity", "Social Service"
),wordwrap, len=36),cex.axis=1,mgp=c(2,1,0),las=1)
#mtext("Treatment",side=1,line=6.5,cex=1.5)
abline(v=0,lty=2,col="red")
for(i in 1:8) segments(tab.approvalexp2[i,2],i,tab.approvalexp2[i,3],i,lwd=1)
segments(tab.approvalexp2[9,2],9,tab.approvalexp2[9,3],9,lwd=1.5)
dev.off()



## load nationwide survey data
nsurvey <- read.dta13("data/nationwidesurvey.dta")
nsurvey$female = 1- nsurvey$gender

## standardize variables
nsurvey_std <- nsurvey

nsurvey_std$scsapproval <- scale(nsurvey_std$credit_systems_approval)
nsurvey_std$distrust <- scale(nsurvey_std$are_chinese_distrusting)
nsurvey_std$avoidfriend <- scale(nsurvey$lower_friends)
nsurvey_std$infotvnewspaper <- scale(nsurvey$infotvnews)
nsurvey_std$infotvnewscmc <- scale(nsurvey$credit_information_TV)
nsurvey_std$femalestd <- scale(nsurvey$female)
nsurvey_std$agestd <- scale(nsurvey$age1)
nsurvey_std$incomestd <- scale(nsurvey$income)
nsurvey_std$ccpstd <- scale(nsurvey$ccp_memberdm)
nsurvey_std$pubemploystd <- scale(nsurvey$pubemploy)
nsurvey_std$pilotstd <- scale(nsurvey$living_in_pilot_city)
nsurvey_std$influencestd <- scale(nsurvey$credit_decision_influence)
nsurvey_std$gtruststd <- scale(nsurvey$government_confidence)
nsurvey_std$fairstd <- scale(nsurvey$credit_score_fairness)
nsurvey_std$edustd <- scale(nsurvey$education)
nsurvey_std$cityruralstd <- scale(nsurvey$city_rural)
nsurvey_std$privacyusestd <- scale(nsurvey$privacyuse)

## main model
mod.approval2 <- lm_robust(scsapproval ~ infotvnewspaper + avoidfriend + distrust + 
                             agestd + edustd + femalestd + incomestd + pubemploystd + pilotstd +
                             cityruralstd + ccpstd + factor(region2), cluster = region,
                           data = nsurvey_std, se_type = "stata")
tab.approval2 <- cbind(mod.approval2$coefficients[12:2],mod.approval2$conf.low[12:2],mod.approval2$conf.high[12:2])


#----------------------Figure 6. Sources of Support for SCSs: Nationwide Online Survey-------------------
pdf("outputmain/approval2.pdf", width=7.9, height=5) 
par(mar=c(3.9, 15, 1, 1) + 1)
plot(tab.approval2[,1],1:11,xlim=c(-0.335,0.335),ylim=c(0.5, 11.5),main = "",ylab = "",pch="",cex.lab=1,cex.axis=1, yaxt="n",xlab="Coefficients of Standardized Variables with 95% CIs")
points(tab.approval2[1:9,1],1:9, pch=1)
points(tab.approval2[10:11,1],10:11, pch=16, cex=1.5)
axis(2,11:1,labels=sapply(c("SCS Info. Source: TV & Newspaper", "Avoid Friends with Bad Credits","Perceived Distrust","Age", "Education", "Female", "Income", "State Employee", "Pilot City", "Urban", "CCP Member"
),wordwrap, len=36),cex.axis=1,mgp=c(2,1,0),las=1)
#mtext("Treatment",side=1,line=6.5,cex=1.5)
abline(v=0,lty=2,col="red")
for(i in 1:9) segments(tab.approval2[i,2],i,tab.approval2[i,3],i,lwd=1)
for(i in 10:11) segments(tab.approval2[i,2],i,tab.approval2[i,3],i,lwd=1.5)
dev.off()


#----------------------Figure 7. Attitude toward Friends with Bad Credits-------------------
nsurvey$friendsfct <- factor(nsurvey$friends_credit_score_influence, labels = c ("No", "Yes", "Don't Know"))

plotfriend <- ggplot(data=subset(nsurvey, !is.na(friendsfct)), aes(x=friendsfct, y = prop.table(stat(count)), fill = friendsfct, label = scales::percent(round(prop.table(stat(count)),2),accuracy = 1L))) + 
  geom_bar() + scale_fill_grey(start=0.8, end=0.2) + guides(fill=FALSE) + expand_limits(y = 0.45) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), panel.border = element_rect(colour = "black", fill=NA, size=0.3))+
  geom_text(stat = 'count',
            position = position_dodge(.9), 
            vjust = -0.5, 
            size = 3) +
  scale_y_continuous(labels = scales::percent_format(accuracy = 1L)) +
  xlab("Look at a Discredited Friend Differently") + ylab("Percentage")

pdf("outputmain/ggfriend.pdf", 3.5, 3.5)
plotfriend
dev.off()

